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Abstract — A conformal dispersive finite-difference time- 
domain (FDTD) method is developed for the study of one- 
dimensional (1-D) plasmonic waveguides formed by an array 
of periodic infinite-long silver cylinders at optical frequencies. 
The curved surfaces of circular and elliptical inclusions are 
modelled in orthogonal FDTD grid using effective permittivities 
(EPs) and the material frequency dispersion is taken into account 
using an auxiliary differential equation (ADE) method. The 
proposed FDTD method does not introduce numerical instability 
but it requires a fourth-order discretisation procedure. To the 
authors' knowledge, it is the first time that the modelling of 
curved structures using a conformal scheme is combined with 
the dispersive FDTD method. The dispersion diagrams obtained 
using EPs and staircase approximations are compared with those 
from the frequency domain embedding method. It is shown that 
the dispersion diagram can be modified by adding additional ele- 
ments or changing geometry of inclusions. Numerical simulations 
of plasmonic waveguides formed by seven elements show that 
row(s) of silver nanoscale cylinders can guide the propagation of 
light due to the coupling of surface plasmons. 



I. Introduction 

It is well known that photonic crystals (PCs) offer unique 
opportunities to control the flow of light [1]. The basic idea 
is to design periodic dielectric structures that have a bandgap 
for a particular frequency range. Periodic dielectric rods with 
removed one or several rows of elements can be used as 
waveguiding devices when operating at bandgap frequencies. 
A lot of effort has been made to obtain a complete and wider 
bandgap. It has been shown that a triangular lattice of air holes 
in a dielectric background has a complete bandgap for TE 
(transverse electric) mode, while a square lattice of dielectric 
rods in air has a bandgap for TM (transverse magnetic) 
mode [2]. The devices operating in the bandgap frequencies 
are not the only option to guide the flow of light. Another 
waveguiding mechanism is the total internal reflection (TIR) 
in one-dimensional (1-D) periodic dielectric rods [3]. It is 
analysed in [3] that a single row of dielectric rods or air holes 
supports waveguiding modes and therefore can be also used as 
waveguide. In [4], the design of such waveguides consisting 
of several rows of dielectric rods with various spacings is 
proposed. 

Recently, a new method for guiding electromagnetic waves 
in structures whose dimensions are below the diffraction limit 
has been proposed. The structures are termed as 'plasmonic 
waveguides' which have an operation of principle based on 
near-field interactions between closely spaced noble metal 
nanoparticles (spacing <C A) that can be efficiently excited 
at their surface plasmon frequency. The guiding principle 
relies on coupled plasmon modes set up by near-field dipole 



interactions that lead to coherent propagation of energy along 
the array. Analogous structures as waveguides in microwave 
regime include periodic metallic cylinders to support propagat- 
ing waves [5], array of flat dipoles which support guided waves 
[6], and Yagi-Uda antennas [7], [8] etc. However, although 
these structures can be scaled to optical frequencies with 
appropriate material properties, their dimensions are limited 
by the so-called diffraction limit A/(2n). On the other hand, 
plasmonic waveguides employ the localisation of electromag- 
netic fields near metal surfaces to confine and guide light in 
regions much smaller than the free space wavelength and can 
effectively overcome the diffraction limit. Previous analysis of 
plasmonic structures includes the plasmon propagation along 
metal stripes, wires, or grooves in metal [9]-[14], and the 
coupling between plasmons on metal particles in order to 
guide energy [15], [16] etc. Such sub wavelength structures 
can also find their applications e.g. efficient absorbers and 
electrically small receiving antennas at microwave frequencies. 
Recently composite materials containing randomly distributed 
electrically conductive material and non-electrically conduc- 
tive material have been designed [17]. They are noted to ex- 
hibit plasma-like responses at frequencies well below plasma 
frequencies of the bulk material. 

The finite-difference time-domain (FDTD) method [18] 
is seen as the most popular numerical technique especially 
because of its flexibility in handling material dispersion as 
well as arbitrary shaped inclusions. In [19], the optical pulse 
propagation below the diffraction limit is shown using the 
FDTD method. Also with the FDTD method, the waveguide 
formed by several rows of silver nanorods arranged in hexag- 
onal is studied [20]. Despite these examples of applying the 
FDTD method for the plasmonic structures, the accuracy of 
modelling has not been proven yet. When modelling curved 
structures, unless using extremely fine meshes, due to the 
nature of orthogonal and staggered grid of conventional FDTD, 
often modifications need to be applied in order to improve 
the numerical accuracy, such as the treatment of interfaces 
between different materials even for planar structures [21], and 
the improved conformal algorithms using structured meshes 
[22] for curved surfaces. 

In addition to the modifications at material interfaces, the 
material frequency dispersion has also to be taken into account 
in FDTD modelling [23]-[25]. However, modelling dispersive 
materials with curved surfaces still remains to be a challenging 
topic due to the complexity in algorithm as well as the intro- 
duction of numerical instability. An alternative way to solve 
this problem is based on the idea of effective permittivities 
(EPs) [26]-[28] in the underlying Cartesian coordinate system, 
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Fig. 1. Comparison of the filling ratio for E y component in FDTD modelling 
of a circular cylinder using (a) staircase approximations and (b) a conformal 
scheme. The radius of circular cylinder is ten cells. 



and the dispersive FDTD scheme can be therefore modified 
accordingly without affecting the stability of algorithm. In this 
paper, we first propose a novel conformal dispersive FDTD 
algorithm combining the EPs together with an auxiliary differ- 
ential equation (ADE) method [18], then apply the developed 
method to the modelling of plasmonic waveguides formed by 
an array of circular or elliptical shaped silver cylinders at 
optical frequencies. Numerical FDTD simulation results are 
verified by a frequency domain embedding method [29]. To 
the authors' knowledge, it is the first time that a conformal 
scheme is combined with the dispersive FDTD method for 
the modelling of nano-plasmonic waveguides. 

II. Conformal Dispersive FDTD Method Using 
Effective Permittivities 

Conventionally, staircase approximations are often used to 
model curved electromagnetic structures in an orthogonal 
FDTD domain. Figure [Ha) shows an example layout of an 
infinite-long cylinder in the free space represented in a two- 
dimensional (2-D) orthogonal FDTD domain. The approxi- 
mated shape introduces spurious numerical resonant modes 
which do not exist in actual structures. On the other hand, 
using the concept of filling ratio, which is defined as the ratio 
of the area of material S2 to the area of a particular FDTD cell, 
the curvature can be properly represented in FDTD domain 
as shown in Fig. Htb), where different levels of darkness 
indicate different filling ratios of material e^. The accuracy 
of modelling can be significantly improved compared with 
staircase approximations, as will be shown in a later section. 

According to [28], the EP in a general form is given by 



£ e ff = £ii(l - n 2 ) +e±n 2 , 



(1) 



where n is the projection of the unit normal vector n along 
the field component as shown in Fig. [2 en and e± are par- 
allel and perpendicular permittivities to the material interface, 
respectively and defined as 



£|| = /£2 + (l-/)ei, 

e± = \f/e2 + (l-f)/ei]-\ 



(2) 
(3) 



where / is the filling ratio of material 82 in a certain FDTD 
cell. 

In this paper we consider the inclusions as silver cylinders, 
which at optical frequencies can be modelled using the Drude 




Fig. 2. Layout of a quarter circular inclusion in orthogonal FDTD grid for 
E y component. The radius of circular cylinder is three cells. 



dispersion model 

e 2 (v) 



so 1 



(4) 



where u p is the plasma frequency and 7 is the collision 
frequency. At the frequencies below the plasma frequency, 
the real part of permittivity is negative. In this paper, we 
assume that the silver cylinders are embedded in the free space 

(£1 e ). 

In order to take into account the frequency dispersion of the 
material, the electric flux density D is introduced into standard 
FDTD updating equations. At each time step, D is updated 
directly from H and E can be calculated from D through the 
following steps. Substitute (0) and ® into (Q]) and using the 
expressions for e\ and £2 ©, the constitutive relation in the 
frequency domain reads 



{to 4 - 2 1 jco 3 - [ 7 2 + (1 - J» 2 + 7 (1 
= [w 4 - 27JW 3 - (7 2 + ujp)u> 2 + -yujpjui 
+/(l-/)(l-n>l]e E. 



(5) 



Using the inverse Fourier transformation i.e. jui — > d/dt, 
obtain the constitutive relation in the time domain as 
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The FDTD simulation domain is represented by an equally 
spaced three-dimensional (3-D) grid with periods Ax, Ay 
and Az along x-, y- and ^-directions, respectively. The time 
step is At. For discretisation of ©, we use the central 
finite difference operators in time (5t) and the central average 
operator with respect to time (/x t ): 

dt A (At) 4 ' dt 3 (Aty 
d d t 



dt 2 



(At) 



dt At 



(7) 
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where the operators S t and ji t are defined as in [30]: 
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Here F represents field components and m x ,m y ,m z are 
indices corresponding to a certain discretisation point in the 
FDTD domain. The discretised Eq. © reads 
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Note that in the above equations we have kept all terms to 
be the fourth-order to guarantee numerical stability. Equation 
( TTOb can be written as 
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The indices m x , m y and m z are omitted from (ITTb since E 
and D are located at same locations. Solving for E n+1 then 
the updating equation for E in FDTD iterations reads 

E n+1 = [6 D n+1 + 6iD n + b 2 B n - 1 + b 3 D n ~ 2 + & 4 D n_3 

- (aiE n + a 2 E n ~ 1 + a 3 E n - 2 + a 4 E n - 3 )]/a , (12) 

with the coefficients given by 
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The computations of H and D are performed using Yee's 
standard updating equations in the free space. Note that if the 
plasma frequency is equal to zero (u p = 0), then (IT2l) reduces 
to the updating equation in the free space i.e. E = D/sq. 

III. FDTD Calculation of Dispersion Diagram 

Applying the Bloch's periodic boundary conditions (PBCs) 
[31]— [36], FDTD method can be used to model periodic 
structures and calculate their dispersion diagrams [37], [38]. 
For any periodic structures, the field at any time should satisfy 
the Bloch theory, i.e. 

= E(d)e^' ka , H(d + a) = H(d)e^' ka , 



E(d- 



(13) 



where d is the distance vector of any location in the com- 
putation domain, k is the wave vector and a is the lattice 
vector along the direction of periodicity. When updating the 
fields at the boundary of the computation domain using FDTD, 
the required fields outside the computation domain can be 
calculated using known field values inside the domain through 
Eq. (fT3l) . Although instead of using real values in conventional 
FDTD computations, the calculation of dispersion diagrams 
requires complex field values, since only one unit cell is 
modelled, the computation load is not significantly increased. 

First we apply the developed conformal dispersive FDTD 
method to calculate the dispersion diagram for 1-D plasmonic 
waveguides formed by an array of periodic infinite-long (along 
z-direction) circular silver cylinders. As shown in Fig. [3j 
the 2-D simulation domain (x-y) with TE modes (therefore 
only E x , E y and H z are non-zero fields) is truncated using 
Bloch's PBCs in ^-direction and Berenger's perfectly matched 
layers (PMLs) [39] in ^-direction. The Berenger's PMLs have 
excellent performance for absorbing propagating waves [39], 
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Fig. 3. The layout of the 2-D FDTD computation domain for calculating 
dispersion diagram for 1-D periodic structures. The inclusion has a circular 
cross-section with radius r and the period of the 1-D infinite structure is a. 



however for evanescent waves, field shows growing behaviour 
inside PMLs. Since the waves radiated by point or line sources 
consist of both propagating and evanescent components, extra 
space (typically a quarter wavelength at the frequency of 
interest) between PMLs and the circular inclusion is added 
to allow the evanescent waves to decay before reaching the 
PMLs. 

The radius of silver cylinders is r = 2.5 x 10 -8 m and 
the period is a = 7.5 x 10 -8 m. The plasma and collision 
frequencies are co p = 9.39 x 10 15 rad/s and 7 = 3.14 x 10 13 
Hz, respectively in order to closely match the bulk dielectric 
function of silver [40]. The FDTD cell size is Ax = Ay = 
2.5 x 10 -9 m with the time step At = Ax/(^/2c) s (where c 
is the speed of light in the free space) according to the Courant 
stability criterion [18]. Although the stability condition for 
high-order FDTD method is typically more strict than the 
conventional one, since the average operator ji t is applied to 
develop the algorithm, we have not found any instability for 
a complete time period of more than 40,000 time steps used 
in all simulations. 

A wideband magnetic line source is placed at an arbitrary 
location in the free space region of the 2-D simulation domain 
in order to excite all resonant modes of the structure within 
the frequency range of interest (normalised frequency / = 
ua/(27rc) e [0 - 0.5]): 



g(t) = e -(^) 2 <e jujt 



(14) 



where to is the initial time delay, r defines the pulse width 
and ou is the centre frequency of the pulse (/ = 0.25). The 
magnetic fields at one hundred random locations in the free 
space region are recorded during simulations, transformed 
into the frequency domain and combined to extract individual 
resonant mode corresponding to each local maximum. For 
each wave vector, a total number of 40,000 time steps are 
used in our simulations to obtain enough accurate frequency 
domain results. 

In order to demonstrate the advantage of EPs and validate 
the proposed conformal dispersive FDTD method, we have 
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Fig. 4. Comparison of the first resonant frequency (transverse mode) at 
wave vector k x = tt/cl calculated using the FDTD method with staircase 
approximations, the FDTD method with EPs and the frequency domain 
embedding method. 



also performed simulations using staircase approximations for 
the circular cylinder, as shown in Fig.[TJa). Figured shows the 
comparison of the first resonant frequency (transverse mode) at 
wave vector k x = n/a of the plasmonic waveguide calculated 
using the FDTD method with staircase approximations, the 
FDTD method with EPs and the frequency domain embed- 
ding method [41]. With the same FDTD spatial resolutions, 
the model using EP shows excellent agreement with the 
results from the frequency domain embedding method, on 
the contrary, the staircase approximation not only leads to 
a shift of the main resonant frequency, but also introduces 
a spurious numerical resonant mode which does not exist in 
actual structures. The same effect has also been found for non- 
dispersive dielectric cylinders [42]. It is also shown in Fig. H] 
that although one may correct the main resonant frequency 
using finer meshes, the spurious resonant mode still remains. 

The problem of frequency shift and spurious modes become 
severer when calculating higher guided modes near the 'flat 
band' region (i.e the region where waves travel at a very low 
phase velocity). Even with a refined spatial resolution, the 
staircase approximation fails to provide correct results (not 
shown). On the other hand, using the proposed conformal 
dispersive FDTD scheme, all resonant modes are correctly 
captured in FDTD simulations as demonstrated by the com- 
parison with the embedding method as shown in Fig. 

According to previous analysis using the frequency embed- 
ding method, the fundamental mode in the modelled plasmonic 
waveguide is transverse mode and the second guided mode is 
longitudinal [41], which is also shown by the distribution of 
electric field intensities in Fig. [6] from our FDTD simulations. 
The higher guided modes are referred to as plasmon modes. 
For demonstration of field symmetries and due to the TE mode 
considered in our simulations, we have plotted the distributions 
of magnetic field corresponding to different resonant modes 
at wave number k x = ir/a as marked in Fig. as shown 
in Fig. [71 Sinusoidal sources for excitation of certain single 
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Fig. 5. Comparison of dispersion diagrams for an array of infinite-long (along 
z -direction) circular silver cylinders calculated using the FDTD method with 
EPs and the frequency domain embedding method. 
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Fig. 6. Normalised total electric field intensities corresponding to (a) 
transverse and (b) longitudinal modes [41] at wave number k x = tr/a as 
marked in Fig. \5\ The structure is infinite along x-direction. 
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Fig. 7. Normalised distributions of magnetic field corresponding to different 
resonant modes at wave number k x = tt/cl as marked in Fig. \5\ (a), (c), 
(d): even modes, and (b), (e): odd modes. The structure is infinite along x- 
direction. (Note that the coordinate has been rotated 90 degrees anti-clockwise 
from Fig. [3] for better presentation of the figure.) 
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mode are used and the sources are placed at different locations 
corresponding to different symmetries of field patterns. All 
field patterns are plotted after the steady state is reached in 
simulations. The modes (a), (c) and (d) in Fig. [7] are even 
modes (relative to the direction of periodicity of the waveguide 
i.e. x-axis), and (b) and (e) are considered as odd modes. 

The above comparison of the simulation results calculated 
using the conformal dispersive FDTD method and the embed- 
ding method clearly demonstrates the effectiveness of applying 
the EPs in FDTD modelling. Furthermore, in contrast to the 
embedding method, the main advantage of the FDTD method 
is that arbitrary shaped geometries can be easily modelled. We 
have applied the conformal dispersive FDTD method to study 
the effect of different inclusions on the dispersion diagrams 
of 1-D plasmonic waveguides. The geometries considered are 
two rows of periodic infinite-long (along z-direction) circular 
silver cylinders arranged in square lattice and a single row 
of periodic infinite-long (along z -direction) elliptical silver 
cylinders. The elliptical element has a ratio of semimajor-to- 
semiminor axis 2:1, where the semiminor axis is equal to the 
radius of the circular element (25.0 nm). For the two rows 
of circular nanorods, the spacing between two rows (centre- 
to-centre distance) is 75.0 nm. The dispersion diagrams for 
these structures are plotted in Figs. [8] and [TO) Comparing the 
dispersion diagrams for a single circular element in Fig. [5] and 
two circular elements in Fig. [8] we can see that the dispersion 



Fig. 8. Dispersion diagram for two rows of periodic infinite-long (along z- 
direction) circular silver cylinders arranged in square lattice calculated from 
conformal dispersive FDTD simulations. 

diagram has been modified due to the change of inclusion. The 
strong coupling between two elements introduces additional 
guided modes to appear in dispersion diagram. Such phe- 
nomenon has also been studied for dielectric (non-dispersive) 
nanorods previously [4]. The distributions of magnetic field 
for selected guided modes as marked in Fig. [8] are plotted in 
Fig. [3 The modes (a), (c) and (d) are even modes while (b) 
and (e) are odd modes. 

The dispersion diagram for a single elliptical element as 
inclusion is shown in Fig. [TO) It can be seen that more guided 
modes appear which is caused by the change of inclusion's 
geometrical shape from circular to elliptical. The frequency 
corresponding to the lowest mode has been lowered due to the 
increase of inclusion's volume. The distributions of magnetic 
field for selected guided modes are plotted in Fig. [TT] The 
modes (a) and (d) are even modes and (b), (c) and (e) are odd 
modes, respectively. 

IV. Wave Propagation in Plasmonic Waveguides 

FORMED BY FINITE NUMBER OF ELEMENTS 

In order to study wave propagations in plasmonic waveg- 
uides formed by a finite number of silver nanorods, we have 
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Fig. 9. Normalised distributions of magnetic field corresponding to different 
guided modes as marked in Fig. [8] (a), (c), (d): even modes, and (b), (e): odd 
modes. The structure is infinite along x-direction. (Note that the coordinate 
has been rotated 90 degrees anti-clockwise from Fig.[3]for better presentation 
of the figure.) 
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Fig. 11. Normalised distributions of magnetic field corresponding to different 
guided modes as marked in Fig.[To] (a), (d): even modes, and (b), (c), (e): odd 
modes. The structure is infinite along x-direction. (Note that the coordinate 
has been rotated 90 degrees anti-clockwise from Fig. [3] for better presentation 
of the figure.) 
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Fig. 10. Dispersion diagram for a single row of periodic infinite-long (along 
z -direction) elliptical silver cylinders calculated from conformal dispersive 
FDTD simulations. 



replaced PBCs in x-direction with PMLs and added additional 
cells for the free space region to the simulation domain. 
The number of nanorods under study is seven. The spacing 
(pseudo-period) between adjacent elements remains the same 
as for infinite structures considered in the previous section. For 
a single mode excitation, we choose the frequency of corre- 
sponding mode from dispersion diagram, and excite sinusoidal 
sources at different locations with respect to the symmetry of 
different guided modes at one end of the waveguides. 

For the plasmonic waveguides formed by different types 
of inclusions, we have chosen certain eigen modes: mode 
Fig - Eta) for a single row of circular cylinders, mode Fig. Od) 
for two rows of circular cylinders, and mode Fig. [Til e) for a 
single row of elliptical cylinders. The distributions of magnetic 
field intensity for different waveguides operating in these 
guided modes are plotted in Fig. [121 The field plots are taken 
after the steady state is reached in simulations. It is clearly seen 
that single guided modes are coupled into these waveguides 
but the excitation of certain modes highly depends on the 
symmetry of field patterns. The energy that can be coupled 
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Fig. 12. Normalised distributions of magnetic field intensity corresponding 
to different guided modes for seven-element plasmonic waveguides formed by 
(a) a single row of circular nanorods (the corresponding eigen mode is shown 
in Fig. Eta)) (b) two rows of circular nanorods arranged in square lattice (the 
corresponding eigen mode is shown in Fig. 03c)) (c) a single row of elliptical 
nanorods (the corresponding eigen mode is shown in Fig. [TTJe)). (Note that 
the coordinate has been rotated 90 degrees anti-clockwise from Fig. [3] for 
better presentation of the figure.) 



into the waveguides also depends on the matching between 
source and the plasmonic waveguide. 

V. Conclusions 

In conclusion, we have developed a conformal dispersive 
FDTD method for the modelling of plasmonic waveguides 
formed by an array of periodic infinite-long silver cylinders 
at optical frequencies. The conformal scheme is based on ef- 
fective permittivities and its main advantage is that since con- 
ventional orthogonal FDTD grid is maintained in simulations, 
no numerical instability is introduced. The material frequency 
dispersion is taken into account using an auxiliary differential 
equation method. The comparison of dispersion diagrams 
for one-dimensional periodic silver cylinders calculated using 
the conformal dispersive FDTD method, the conventional 
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dispersive FDTD method with staircase approximations and 
the frequency domain embedding method demonstrates the 
accuracy of the proposed method. It is shown that by adding 
additional element or changing the geometry of inclusions, the 
corresponding dispersion diagram can be modified. Numerical 
simulations of plasmonic waveguides formed by seven ele- 
ments show that the eigen modes in infinite structures can be 
excited but highly depend on the symmetry of field patterns 
of certain modes. Further work includes the investigation 
of the effects of different number of elements in plasmonic 
waveguides on guided modes, and the calculation of group 
velocity of different modes propagating in these waveguides. 
Although results presented in this paper have been focused 
at optical frequencies, with future advances in microwave 
plasmonic materials, novel applications can be found in the 
designs of small antenna and efficient absorbers. 
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